The following article is Open access

A Unified Spectroscopic and Photometric Model to Infer Surface Inhomogeneity: Application to Luhman 16B

and

Published 2022 July 12 © 2022. The Author(s). Published by the American Astronomical Society.
, , Citation Michael K. Plummer and Ji Wang 2022 ApJ 933 163 DOI 10.3847/1538-4357/ac75b9

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

0004-637X/933/2/163

Abstract

Extremely large telescopes (ELTs) provide an opportunity to observe surface inhomogeneities for ultracool objects including M dwarfs, brown dwarfs (BDs), and gas giant planets via Doppler imaging and spectrophotometry techniques. These inhomogeneities can be caused by star spots, clouds, and vortices. Star spots and associated stellar flares play a significant role in habitability, either stifling life or catalyzing abiogenesis depending on the emission frequency, magnitude, and orientation. Clouds and vortices may be the source of spectral and photometric variability observed at the L/T transition of BDs and are expected in gas giant exoplanets. We develop a versatile analytical framework to model and infer surface inhomogeneities that can be applied to both spectroscopic and photometric data. This model is validated against a slew of numerical simulations. Using archival spectroscopic and photometric data, we infer starspot parameters (location, size, and contrast) and generate global surface maps for Luhman 16B (an early T dwarf and one of our solar system's nearest neighbors at a distance of ≈2 pc). We confirm previous findings that Luhman 16B's atmosphere is inhomogeneous with time-varying features. In addition, we provide tentative evidence of longer timescale atmospheric structures such as dark equatorial and bright midlatitude to polar spots. These findings are discussed in the context of atmospheric circulation and dynamics for ultracool dwarfs. Our analytical model will be valuable in assessing the feasibility of using ELTs to study surface inhomogeneities of gas giant exoplanets and other ultracool objects.

Export citation and abstract BibTeX RIS

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. Introduction

Imaging ultracool object surface features such as star spots, clouds, and vortices will help us to better understand habitability, climate, and atmospheric dynamics. The advent of 30 m class extremely large telescopes (ELTs) and the James Webb Space Telescope (JWST) present vast opportunities in the fields of ultracool object imaging. ELTs will be unique in providing high signal-to-noise (S/N) and high resolution spectral data while JWST will provide high S/N photometric and low-to-medium resolution (R ≤ 5000) spectral data. These developments will dramatically improve our ability to map the global surfaces of ultracool objects such as M dwarfs, brown dwarfs (BDs), and gas giant exoplanets.

1.1. Star Spots and M Dwarfs

M dwarf stars comprise the most numerous spectral class in the Milky Way. Due to M dwarfs' small size, they are attractive targets for studies attempting to find habitable terrestrial planets because orbiting planets create correspondingly larger and more easily detected signatures with both the transit and radial velocity (RV) methods. Small rocky planets may be over 3 times more common around M dwarfs than FGK stars (Mulders et al. 2015). The overabundance was later confirmed in comparing RV and transiting surveys (Tuomi et al. 2019; Sabotta et al. 2021; Zink et al. 2020).

Despite the prevalence of planets around M stars, obstacles exist for declaring M dwarf stellar systems ideal for planetary habitability. M dwarfs' low temperatures (compared to G dwarf stars like our Sun), mean planets located at distances where liquid water could exist on a planet's surface, approximately 10% the distance between Earth and the Sun, would experience extensive stellar activity including flares (Kasting et al. 1993; Kopparapu et al. 2013). Fast-rotating M dwarfs, and particularly cooler, later and fully convective stars, are far more likely to flare than their hotter counterparts, which possess radiative zones (Günther et al. 2020). Furthermore, flares and even superflares have been observed in ultracool objects with spectral classes as cool as L5, demonstrating that intense magnetic activity continues beyond the M-class range (Paudel et al. 2018, 2020).

Flares and coronal mass ejections associated with star spots could potentially bathe nearby, likely tidally locked, planets in X-ray and ultraviolet (UV) radiation and charged particles, stripping the planets' nascent atmospheres. However, it remains uncertain whether this radiation will inhibit or catalyze the emergence and growth of burgeoning life (Scalo et al. 2007; Tarter et al. 2007; Zendejas et al. 2010). Laboratory experiments suggest that stellar flares could provide the UV light necessary to photochemically produce the prerequisite molecules for abiogenesis (Rimmer et al. 2018).

Understanding starspot coverage as a proxy for magnetic activity and stellar flares is necessary to reveal how these energetic events may influence the long-term habitability of planetary systems. Late M dwarfs have been observed to have star spots preferentially located within their polar latitudes which differs from the predominantly equatorial-forming Sun spots (Strassmeier & Rice 1998; Järvinen et al. 2008; Barnes et al. 2015, 2017). TESS observations have also provided preliminary evidence of white-light flares preferentially occurring at polar latitudes (between 55° and 81°) in fully convective M dwarf stars (Ilin et al. 2021). Future M dwarf imaging will be critical to confirming theories on starspot formation and migration, and will also be necessary to understand and quantify habitability in these star systems.

1.2. BDs at the L/T Transition

Lower in mass than M dwarfs, BDs are substellar objects below the hydrogen-burning minimum mass (M ≲ 70 MJ ) (Chabrier et al. 2000; Dupuy & Liu 2017). They gradually cool over billions of years from early post-formation temperatures consistent with late M dwarfs, through L, T, and Y substellar classes. Across this span, the spectral features transition quite dramatically in JK color from red, cloudy L Dwarfs to relatively blue, cloud-free mid-T Dwarfs due to potentially two primary mechanisms related to the decreasing temperature: (1) Clouds composed of refractory materials such as iron (Fe), corundum (Al2O3), enstatite (MgSiO3), and forsterite (Mg2SiO4) condense and precipitate out over the L/T transition leading to these materials seemingly being sequestered below the photosphere (2) carbon monoxide (CO) is reduced to methane (CH4) resulting in J-band brightening and H and K suppression due to the ensuing CH4 absorption (Tsuji et al. 1996; Jones & Tsuji 1997; Noll et al. 2000; Allard et al. 2001; Marley et al. 2002; Burrows et al. 2006; Cushing et al. 2008).

BDs at the L/T transition appear to demonstrate cloud patchiness, which may explain the source of variability. It was initially theorized that the condensation of the metal and silicate clouds along with robust tropospheric convection leads to gaps in cloud coverage (Ackerman & Marley 2001; Burgasser et al. 2002; Marley et al. 2010). More recently, it has been argued that a more likely scenario is that condensation results in patchy regions of varying cloud height and optical thickness (Radigan et al. 2012; Apai et al. 2013; Buenzli et al. 2014). Inhomogeneous atmospheres paired with the fact most BDs spin rapidly with average rotational speeds exceeding 20 km s−1 for mid-L dwarfs, naturally leads to the conclusion that ultracool objects at the L/T transition will experience large photometric and spectroscopic variability (Reiners & Basri 2008).

Previous observations indicated that, although objects at the L/T transition were variable, they were perhaps no more variable than other BDs (Enoch et al. 2003; Buenzli et al. 2014; Metchev et al. 2015). However, Radigan et al. (2014) deduced that ${39}_{-14}^{+16} \% $ of L/T transition stars were highly variable (>2%) while also demonstrating strong variability to be rare for BDs outside this range. Eriksson et al. (2019) determined a similar value of ${40}_{-19}^{+32} \% $ for the proportion of strongly variable L/T transition BDs.

Future opportunities exist to shed more light on the exact nature of the L/T transition and BD atmospheric dynamics. ELTs and the JWST will provide high quality data that could potentially be used to map the surfaces of numerous nearby BDs. These maps very well may reveal the existence of atmospheric structure including zonal jets and vortices as well as the true source of L/T transition variability. With this information, we can further constrain atmospheric formation, structure, and dynamics.

Not only will we be able to map these ultracool objects' surface temperatures, we may also be able to map chemical tracers via Doppler imaging. This could allow us detect regions of chemical disequilibrium (Section 5.3) originating from vertical atmospheric flow. Studying such behavior across spectral classes can allow us to fully understand the development of thermal and hydrodynamic structures with ultracool objects.

1.3. Clouds in Gas Giant Exoplanet Atmospheres

Directly imaged gas giant exoplanets typically form at large distances from their host stars and are weakly illuminated, leading to targets similar to BDs in terms of observational opportunities. Like their more massive BD cousins, gas giant exoplanets possess clouds (Madhusudhan et al. 2011; Currie et al. 2011, 2014; Marley et al. 2012; Skemer et al. 2014) and chemical disequilibrium (Barman et al. 2011, 2015; Konopacky et al. 2013; Skemer et al. 2014; Currie et al. 2014; Lavie et al. 2017; Mollière et al. 2020; Wang et al. 2020).

High S/N and high spectral resolution data from ELTs will be essential to creating gas giant surface maps via Doppler imaging in the future. These maps may have the ability to determine the extent of cloud patchiness as well as the occurrence of planetary-scale structures such as zonal jet bands, vortices, and persistent storms. As is later discussed in Section 5.3, Doppler imaging may also be used to identify regions of convective activity through chemical disequilibrium tracers.

1.4. Doppler Imaging

Doppler imagining is an important spectroscopic method for studying surface inhomogeneity. The history of Doppler imaging goes back almost 70 yr when Deutsch (1958) first used spectral line profile (LP) deviations to map the magnetic fields and abundance anomalies of Ap stars; this method was further formalized in Deutsch (1970) and Falk & Wehlau (1974). Ap stars continued to be targeted for surveys mapping their spots and chemical compositions due the class possessing strong surface inhomogeneities with century-long lifetimes (Khokhlova 1976; Goncharskii et al. 1977, 1982). Following this work, Vogt & Penrod (1983) conducted a groundbreaking study involving the Doppler imaging of the spotted star, HR1099.

These early Doppler imaging efforts suffered from the fact that stellar surface imaging is an ill-posed problem, the results were often degenerate with more than one solution for a given set of data (Piskunov et al. 1990). Solutions to this problem arose in applying maximum entropy methods (Vogt et al. 1987) and Tikhonov regularization (Piskunov et al. 1990) to infer the smoothest stellar surfaces.

Since these developments, studies have generated Doppler images of giants and sub-giants (Strassmeier et al. 1998, 1999; Donati 1999; Strassmeier 1999), pre-main-sequence stars (Hatzes 1995; Strassmeier & Rice 1998; Finociety et al. 2021), K Dwarfs (Collier Cameron & Unruh 1994; Barnes et al. 2001, 2004; Zaire et al. 2021), and even M dwarfs (Barnes & Collier Cameron 2001; Barnes et al. 2004, 2017). Notably, Crossfield et al. (2014) conducted the first Doppler imaging of a BD on the nearby binary system Luhman 16; these results were reprocessed and replicated in Luger et al. (2021).

1.5. Photometry and Low-resolution Techniques

Photometric and low-resolution spectroscopic (spectrophotometric) campaigns have led to improved understanding of BD and directly imaged gas giant exoplanet atmospheres and clouds. Atmospheric structure can be investigated by monitoring ultracool object rotational modulation (Buenzli et al. 2012, 2015a, 2015b; Apai et al. 2013; Karalidi et al. 2016; Yang et al. 2016; Schlawin et al. 2017; Biller et al. 2018; Zhou et al. 2020a; Manjavacas et al. 2021). The Hubble Space Telescope (HST) Cloud Atlas program has focused on using these methods and performed extensive spectrophotometric studies of BDs, gas giants, and planetary-mass objects (Lew et al. 2016, 2020b, 2020a; Manjavacas et al. 2018, 2019b, 2019a; Zhou et al. 2018, 2019, 2020b; Miles-Páez et al. 2019). Moving forward, these methods can potentially be integrated with techniques such as Doppler imaging to create more complete surface maps for ultracool objects. The unified model presented in this paper seeks to address this need by combining photometric techniques with Doppler imaging data sets.

1.6. Outline

In this paper, we develop analytical and numerical frameworks for Doppler imaging ultracool objects using the Euler−Rodrigues formula (Section 2). We internally validate our analytical method versus a numerical stellar/substellar surface model (Section 3). Our model is then externally validated against spectroscopic and photometric archival data from Crossfield et al. (2014) and Buenzli et al. (2015b) (Section 4). Based on our results, we discuss the implications for Doppler imaging of ultracool objects and future avenues of research (Section 5).

2. Methods

2.1. Least Squares Deconvolution

LPs have a higher S/N than individual spectral lines and are therefore more effective in identifying surface inhomogeneities. One common method of computing LPs for spectra uses the cross-correlation function; however, for fast-rotating objects, least squares deconvolution (LSD) can be a mathematically rigorous technique. LSD is an attractive option for both spectroscopic and spectropolarimetric applications because it can create high-precision, averaged LPs, ideal for Bayesian inference and Doppler imaging (Kochukhov et al. 2010).

To understand how LSD works, we must consider how the target's observed spectral lines are broadened through multiple processes including rotation and instrumentation. We can account for these processes mathematically by convolving the source spectrum S with a broadening and instrument kernel, represented by Z. This results in the observed spectrum Y0 determined as such,

Equation (1)

where ∗ is the convolution operator. In Equation (1), the convolution operator broadens the source spectrum, S, by the width of the broadening kernel, Z, creating a mathematical representation of the observed spectrum.

In most cases, we do not have perfect information and cannot directly observe the broadening kernel. The source spectrum is only estimated based on prior observations and models. The parameters we wish to determine are often innate to the source spectrum, but our primary observable is the observed spectrum, Y0. We estimate the source spectrum with a template spectrum, F, which can be computed from various models based on known or estimated parameters such as effective temperature (Teff), log(g), and metallicity.

We conducted sensitivity tests for both effective temperature and log(g) with a fixed template spectra (Teff = 1500 K, log(g) = 5.0, [Fe/ H] = 0), in the range used in Crossfield et al. (2014). We created synthetic observed spectra with varying effective temperature (±50 K) and log(g) (±0.5); available spectra in this temperature regime did not have multiple metallicities available. This synthetic spectra was implemented, along with the fixed template, into our numerical model (see Section 2.3) to create LPs with a spot of known latitude, longitude, contrast, and radius. We then retrieved these parameters via Bayesian inference. A temperature difference of 50 K resulted in the following deviations in the retrieved parameters: latitude (1.24σ), longitude (0.57σ), radius (0.57σ), and contrast (0.19σ). A log(g) difference of 0.5 resulted in the following: latitude (1.68σ), longitude (0.096σ), radius (0.9σ), and contrast (0.16σ). As can be seen, longitude, radius, and especially contrast were fairly robust to temperature and gravity deviations while latitude was less robust.

We can replace S with F in Equation (1) as follows:

Equation (2)

In this equation, Y0 is comprised of m elements with spacing dependent on the spectral sampling of the observational instrument. In practice, solving Equation (2) requires expanding F into a line mask, M, via an m × n Toeplitz matrix as described in Donati et al. (1997). This provides us with the following form:

Equation (3)

Similar to Crossfield et al. (2014), we use LSD to convert the observed spectrum into an LP in RV space. These LPs can then be compared to either a mean LP from the observation or to a modeled LP to infer surface inhomogeneities. The specific LSD method used in this paper was adopted from Wang et al. (2017, 2018) and Pai Asnodkar et al. (2022), and is best represented by the following:

Equation (4)

where Z is the computed LP comprised of n data points and S is an m × m diagonal error matrix with Sii elements equal to 1/σii where σii is the measurement error at each epoch. Λ is a scalar regularization parameter while R is the associated regularization Toeplitz matrix, which is outlined in Donatelli & Reichel (2014).

2.2. Analytical Model

2.2.1. Scheme

Doppler tomography examines the changes in a target's LP during a transit of its viewable hemisphere. Our analytical method for modeling stellar and substellar surfaces derives inspiration from Doppler tomographic studies of eclipsing binary systems (Albrecht et al. 2007; Wang et al. 2018) and transiting exoplanets (Hirano et al. 2011; Johnson et al. 2014). Spectral deviations created by a transiting companion or surface inhomogeneity blocking a portion of the target's visible surface can be measured using a high-S/N technique such as LSD.

These spectral deviations create Doppler shadows when they are plotted in RV space. With this technique, companion stars and planets will create linear shadows when they transit the target. In this article, we make the unique contribution of considering that surface inhomogeneities create curved Doppler shadows, if the target is inclined with respect to the viewer, as can be seen in Figure 1. For the hemisphere inclined toward the viewer, higher latitudes spots are more extensive in their viewable phase.

Figure 1.

Figure 1. Surface inhomogeneity trajectories and corresponding LP deviation plots for example stellar/substellar object with υ = 50 km s−1 and i = 60°. (Left) Equatorial Spot (l = 0°). (Middle) Midlatitude spot (l = 40°). (Right) Near-polar spot (l = 75°). It should be noted that a transiting planet would create a linear signature vs. the curved signatures seen in the figure.

Standard image High-resolution image

We represent the unperturbed stellar/substellar disk as a rotationally broadened LP. Surface inhomogeneities are modeled as Gaussian distributions whose width depends on radius and viewing angle. These Gaussian distributions are either subtracted from (dark spot) or added to (bright spot) the rotationally broadened LP to create an analytical model for each epoch. An example of this process can be seen in Figure 2 with dark spots subtracting from the rotationally broadened LP.

Figure 2.

Figure 2. LPs of the spotted target with an inclination, i = 60°, and viewing longitudes, α = −45°,0°, +45°, with an equatorial latitude, radius of 10°, and contrast of 0.5. The black line denotes the unperturbed rotationally broadened LP. (Inset) Corresponding demonstration of a spot at three different longitudes on stellar disk.

Standard image High-resolution image

2.2.2. Unperturbed LP

We approximate the unperturbed LP with a rotational broadening kernel. Rotation is the dominant spectral line broadening mechanism for fast-rotating astrophysical objects ($\upsilon \ \sin \ i\gt 20\ \mathrm{km}\,{{\rm{s}}}^{-1}$, where υ is the equatorial rotation speed of the object and i is its inclination).

A suitable expression for this rotational broadening kernel is described in Gray (2008) and adapted by Pai Asnodkar et al. (2022) as follows:

Equation (5)

G(υi , t) is the rotational broadening kernel applied over velocities υi and epochs, t, while υ is the RV of the star with respect to the solar system. The constants c1 and c2 are computed as follows:

Equation (6)

and

Equation (7)

where epsilon is the linear limb darkening coefficient describing how the center of the stellar or substellar target appears brighter than the edge.

As the linear limb darkening law represents a poor approximation of the limb darkening phenomenon for stellar and substellar targets, we conducted a comparison between the rotational broadening kernel in Equation (5) and one created via higher-order Claret limb darkening. Using Scipy's least squares tool (Virtanen et al. 2020), we fit Equation (5) to an LP created with Claret's limb darkening coefficients for a 1500 K, log(g) = 5.0, solar metallicity, substellar object (Claret et al. 2012). Both profiles were also convolved with an instrument profile. The resulting profiles agreed within 0.2%, demonstrating the suitability of Equation (5) for representing ultracool stellar and substellar disks.

2.2.3. Gaussian Spots

Surface inhomogeneities can be described by Gaussian distributions centered on an RV value that is a function of the longitudinal viewing angle with respect to the observer. The Gaussian distribution's width is set by the surface inhomogeneity's longitudinal radius. The magnitude of the Gaussian distribution is a product of its area and contrast.

As the target star or substellar object rotates, the Gaussian distribution will transit across the LP. The range of the transit in RV space depends on the inhomogeneity's latitude. As demonstrated in Figure 1, an equatorial inhomogeneity will transit the entire rotationally broadened LP while an inhomogeneity at mid-to-polar latitudes will transit an RV band commensurate with its latitude. Conversely, a polar spot on an inclined target may not rotate out of view such that while its RV domain is limited by its latitude, its signature will be viewable throughout the entire phase space.

We follow the examples of Crossfield et al. (2014), Karalidi et al. (2015), and Karalidi et al. (2016), by modeling surface inhomogeneities as circles or ellipses. These inhomogeneities are initially modeled as two-dimensional Gaussian ellipses with the longitudinal and latitudinal widths adjusted for viewing angle by Lambert's cosine law,

Equation (8)

which states that the observed intensity of an Lambertian object decays according to a cosine function where θ is the angle of the observer with respect to the normal. Assuming a circular Gaussian spot, the longitudinal (rL ) and latitudinal (rl ) radii then become

Equation (9)

where r is the radius of the Gaussian spot with units of degrees and α and β are the longitudinal and latitudinal viewing angles. These viewing angles are measured from the sub-observer centroid of target to the center of the Gaussian spot. They are computed by first assuming that the target is a sphere with radius equal to unity and the following spherical coordinates:

Equation (10)

The spot's latitude and longitude are then converted into Euler coordinates, $\vec{x}=(x,y,z)$,

Equation (11)

We define the rotational angle (ξ) as a function of the inclination (ξ = 90°–i). The transformation matrix (T) is composed of rotation parameters (a, b, c, d). Computing these rotation parameters, we encounter the components of the unit vector $\vec{k}=({k}_{x},{k}_{y},{k}_{z})$, but as we will only rotate along the y-axis by varying the inclination, kx = kz = 0 (Shuster 1993).

Equation (12)

The rotated Euler coordinates $(\vec{x^{\prime} })$ can now be computed with the following transformation (Shuster 1993):

Equation (13)

where T is defined as follows (Shuster 1993):

Equation (14)

The Euler−Rodrigues formula is used to account for the coordinate transformation necessitated by any orbital inclination that is not 90°. The transformation results in new coordinates ($\vec{x^{\prime} }$), which can then be converted back into the longitudinal and latitudinal viewing angles, α and β. Using Equation (9), the longitudinal radius, rL , is then converted into RV space for use in the LP; this will be the standard deviation (σ) of the Gaussian spot.

Using the Euler−Rodrigues formula to compute a spot's position with respect to the observer's viewing angle creates a set of coordinates which, depending on the duration of observations, could include spot positions encompassing an entire revolution. In subsequent analysis, we will compute the longitudinal and latitudinal viewing angles for one complete revolution. Depending on the target's inclination, along with the spot's latitude, the spot may have instances where it is not viewable for periods of time. In the opposite case, for positive inclinations (0° < i < 90°), if the latitude is greater than inclination (l > i), spots in the Northern Hemisphere (hemisphere tilted toward the observer) will remain in view for the entire revolution as they circle the North Pole. Each of these cases must be taken into account and spot-induced LP deviations will be set to zero when the spot is not within view of the observer.

The magnitude of the Gaussian spot is a function of the spot's area and contrast. The area (aspot) is computed as a fraction of the star's unit disk:

Equation (15)

This leads us to the Gaussian distribution for the surface inhomogeneity:

Equation (16)

where

Equation (17)

and

Equation (18)

RVspot is the spot centroid RV for each epoch and RV ranges from −3σ to 3σ to capture the continuum well beyond the expected rotationally broadened signal at ±1σ for normalization purposes. The term v sin i/90° is used to convert from longitude to velocity space.

Contrast varies in the range [−1, 1] with −1 being twice as bright as the average surface, 0 demarcating the average surface temperature, and 1 denoting a perfectly dark spot with 0 K temperature and no photon emission. It can be seen that a negative contrast will lead to more flux added to an LP centered at the spot's location while positive contrast will lead to a subtraction of flux centered at the spot's location. The spot's maximum achievable magnitude is simply the product of its area and its contrast.

For each epoch, the computed Gaussian spot in the RV domain is then added or subtracted to the rotational broadening LP. Figure 2 demonstrates the analytical method with a starspot moving across a 60° inclined disk with an equatorial latitude; −45°,0°, and +45° longitudes; 10° radius; and 0.5 contrast.

2.2.4. Analytical Model Photometry

Although the analytical model was initially created to process observed spectra into LPs, it can be adapted to create photometric light curves integrating the total flux at each phase of the target's rotation. The LSD method incorporates a broad spectral band, so we can assume a featureless and normalized rotational broadening kernel has a summed magnitude approximately the equal to unity at each epoch. Spots with positive contrast (dark spots) will subtract from the total flux at each epoch while spots with negative contrast (bright spots) will add to the total flux.

The produced light curves exhibit a rotationally broadened shape with a sharp cutoff as a spot enters and exits the stellar limb. This sharp cutoff leads to issues with inferring spot parameters (particularly when applied to real observations such as Luhman 16B in Section 4.5). To fix this issue, when the analytical code is applied to real-world data, a 1D Gaussian filter (from the SciPy library) is applied to both analytical and numerical photometric light curves (Virtanen et al. 2020).

2.3. Numerical Model

While an analytical starspot model has clear advantages in terms of elegance and computational speed, a numerical model for stellar and substellar photospheric surfaces is necessary to verify the validity of the analytical approach. Building a numerical model for the star also provides a convenient option for creating global maps of the star's surface to verify each spot's location as well as their motion (ensuring that the motion computed by the analytical method is indeed physical and passes the so-called visual check).

To create the numerical model, the star's surface is divided n times longitudinally and n/2 times latitudinally to create n2/2 cells. To speed computation, n longitudinal cells are implemented so that each longitudinal cell directly maps to an RV value. For example, a longitude of +90° directly maps to an RV of + υsin i. Likewise, a longitude of −90° maps to an RV of −υ sin i.

It is also necessary to consider the unseen hemisphere. The observer only has a view of one hemisphere, so the flux at any point located at an angular distance greater than 90° from the sub-observer point (α = 0°, β = 0°), must be set to zero.

All the points within the 90° radius are initially set to 1 before being multiplied by a linear limb darkening law (to align with the linear limb darkening used in the rotational broadening equation in Section 2.2 as well as to match the practice in Crossfield et al. 2014). The insets of Figures 1 and 2 demonstrate the numerical code's ability to produce limb-darkened stellar surfaces.

Star spots can then be added to the numerical model using a bivariate Gaussian distribution with the longitudinal and latitudinal radii (standard deviations) calculated with Equation (9). As in the analytical technique, the Euler−Rodrigues formula is used to transform the spot's coordinates based on the target's inclination for one entire rotation. The spots are then added or subtracted (bright spots are added while dark spots are subtracted) to the limb-darkened image.

The numerical code can be implemented to generate a dynamic LP that accounts for the rotation of spots into and out of view as well as the target's inclination. As mentioned above, the target's longitude directly maps to the RV in the LP. Figure 3 demonstrates the starspot's effect on the LP. Both the numerical and analytical codes lead to a smaller absorption profile in the redshifted hemisphere, which aligns with our expectations that the LP should experience a blueshift when a surface inhomogeneity is in the Eastern Hemisphere (the hemisphere rotating away from the observer).

Figure 3.

Figure 3. (Top) Analytical and numerical LPs with starspot in the positive RV domain indicating a small blueshift. (Bottom) Residuals between the analytical and numerical model LP results with a maximum amplitude of 2.3 × 10−5. The small residual values indicate an excellent match between the two different techniques. (Inset) Surface map of stellar surface with spot at (l = 30°, L = 20°) with a 15° radius and contrast of 0.5 with i = 85°.

Standard image High-resolution image

Similar to the analytical model, the numerical code can also be used to generate photometric light curves. This is accomplished by summing up the total flux across the numerical map at each epoch. This will create a value that far exceeds unity, so normalization is necessary. To accomplish this normalization, a featureless numerical map of the star or substellar object is generated using the selected limb darkening law. The sum of flux at each epoch is then divided by this featureless sum of flux to create a normalized value for each measurement. A light curve closely matching the analytical photometry is produced, as seen in Figure 4.

Figure 4.

Figure 4. Photometric light curve comparison for analytical and numerical methods. Light curves include integrated flux summed across the entire disk for each rotational phase angle. A flux reduction of >3% due to the modeled spot can be seen centered at 180° through the rotational period. Computed residuals between the analytical and numerical models are on the ordered of 10−6. (Inset) Surface map of stellar surface with spot at (l = 0°, L = 180°) with a 20° radius and contrast of 0.3 with i = 85°.

Standard image High-resolution image

A similar issue arises with the numerical curve as that seen with the analytical curve, the spot's curve is not smooth. We apply the same solution and use a 1D Gaussian filter to smooth the curve leading to a more physical appearing light curve that also aids in fitting data with Bayesian inference.

3. Model Validation

We validate the analytical model (Section 2.2) against the numerical mode (Section 2.3). This is necessary to show that the two models are approximately interchangeable, so it follows the analytical model can be used as a replacement for the numerical method when performing potentially computationally expensive posterior inferences.

For the 1-Spot case, the analytical and numerical models' spectroscopic and photometric results match both qualitatively and quantitatively. Figure 3 directly compares LPs generated from the analytical and numerical methods. Residual errors are small on the order of 10−5. A comparison of the photometric light curves created by the analytical and numerical methods can be seen in Figure 4. Again, the residual is small with error on the order of 10−6.

Expanding our test cases, an injection-retrieval technique is used. Two spots are assigned to a stellar surface in the numerical model with each spot possessing a latitude, longitude, radius, and contrast to create a simulated observation and associated LP. Bayesian inference via nested sampling is then used to retrieve the two spots' locations with our analytical method.

We adopt a similar technique to Karalidi et al. (2015) where we explore the angular distance, in terms of latitudinal (Δl) and longitudinal (ΔL) displacement, between the two spots. Both spots must be inferred simultaneously using a 2-Spot model as subsequent retrievals with a 1-Spot model will only infer the spot with the largest LP deviation unless restrictive priors are applied. In terms of priors, we assume uninformative priors for the first spot for both longitude (L ± 180°) and latitude (l ± 90°). For the second spot, we assume either a greater latitude or longitude than the first to overcome the multimodality of the solution. These tests will be accomplished for varying stellar inclinations with respect to the observer and demonstrate improved results compared with those of Karalidi et al. (2015) in terms of the spatial difference required to resolve to stellar/substellar surface features.

3.1. Nested Sampling and Inference

We use the open-source Python module, Dynesty, to infer spots' locations, size, and brightness using our analytical model as well as assumed priors (Speagle 2020). For our purposes, performing Bayesian inference on stellar/substellar surfaces using Markov Chain Monte Carlo (MCMC) methods does not necessarily result in optimal or easy-to-decipher posterior probabilities due to the multimodal nature of star spots and BD atmospheric features once one moves past a 1-Spot model. For this reason, a nested sampling approach was adopted. Dynesty not only functions well with multimodal data, it also computes the marginal likelihood, or in nested sampling terms, the evidence (Speagle 2020).

John Skilling developed nested sampling to solve for the evidence by computing the likelihood and prior mass volume through series of nested volumes of likelihood (Skilling 2004, 2006). Live points are distributed throughout the prior volume and at each iteration the live point with the lowest likelihood is replaced with a new live point with a greater likelihood; the discarded live points (dead points) are then used to compute the evidence and posterior distribution (Skilling 2004, 2006).

3.2. Spectroscopic Validation

In this section, we will perform a sequence of validation tests, each resolving two spots with varying angular separation and inclination.

3.2.1. Test Cases Based on Karalidi et al. (2015)

To accomplish a systematic exploration of the model's robustness, we adopt a similar approach to that of Karalidi et al. (2015), which validated Aeolus, a photometric MCMC-based substellar surface mapping code. Karalidi et al. (2015) tested their code by varying the angular distance in terms of latitude and longitude for two spots with an approximately infinite S/N. Spot coordinates were input and the model was used to retrieve both spots location. For this work, the numerical model was implemented with input spot locations, sizes, and contrasts to generate simulated observed LPs. We then performed a Bayesian inference using Dynesty to determine the spot locations.

As in Karalidi et al. (2015), we assume a larger spot (fixed at the equator at the observer's latitude) with a radius of 9° and contrast of 0.3. The smaller spot (with variable location) has a radius of 5° and contrast of 0.6. For this validation test, $v\ \sin \ i=50\ \mathrm{km}\,{{\rm{s}}}^{-1}$ was selected based on the rotational range recommended by Vogt et al. (1987), and the inclination is nearly edge-on with i = 85°.

For each case, the input locations are retrieved within 1σ for both spots. As can be seen in Figure 5, the analytical model in this work improves on the angular distance that can be resolved: both the smaller angular separations of 24° longitudinally and 20° latitudinally can be resolved.

Figure 5.

Figure 5. Spectral injection-retrieval comparison. (a) Figure from Karalidi et al. (2015). Spots are input and retrieved with longitudinal (ΔL) and latitudinal (Δl) displacements of 24°, 74°, 20°, and 80°, respectively. Separate spots are not retrieved for the ΔL = 24°, Δl = 20°, and Δl = 80° cases. (b) Input and retrieval maps created using numerical and analytical models in this work. For each case, both input spots are retrieved within 1σ. Simulated LPs were created using the numerical model based on spot locations, and these locations were inferred using nested sampling with the analytical model likelihood implementation.

Standard image High-resolution image

We note that the reader may notice in Figure 5 the angular distances from Karalidi et al. (2015) do not appear to exactly match the distances from the analytical model from this work. This is an artifact relating to how the maps in Karalidi et al. (2015) are constructed with two-dimensional circles while the maps generated in this work are effectively three dimensional due to the Euler−Rodrigues transformation. For example, an equatorial spot traveling 10° across a stellar disk will appear to travel less because of the larger radial motion along the line of sight than the tangential motion.

Another important difference between the two models is this paper's adoption of Lambert's cosine law. In Figure 5, the effects of the law can be seen in how the spot's shape is modified near the stellar limb to account for the foreshortening of the received flux for a spot seen at high viewing angles.

An example of the input simulated LP data and the retrieved 2-Spot analytical model for the ΔL = 24° case can be seen in Figure 6. Deviations for both spots can be seen to transit the LP during the early and late phases of the rotation. The 1% noise level can also be observed in the simulated data.

Figure 6.

Figure 6. Example LPs for longitudinally displaced (ΔL = 24°) case. Includes 10 epochs, representing one rotation (represented in the figure by the phase angle). Dashed red lines represent input numerical spectra with a 1% noise level while the solid blue lines denote the 2-Spot analytical model. LPs are offset by 0.02 for comparison, but are normalized to unity within the code. The deviations due to surface inhomogeneities can be seen in both the input and retrieved LPs.

Standard image High-resolution image

3.2.2. Comprehensive Grid Test

Expanding on the comparison shown in Figure 5, a more exhaustive parameter space validation was conducted with longitude and latitude displacements varying between 0° and 60°. A realistic noise level of 1% was selected and radii and contrast were set to 20° and 0.3, respectively, to reflect previous BD Doppler imaging results for cloud system radii from Crossfield et al. (2014) and Luger et al. (2021) as well as typical spot temperatures observed in the Sun's photosphere. Inclinations of 30° and 60° were selected to explore the code's ability to resolve spots and break degeneracies (seen at edge-on inclinations) at different viewing angles. As above, the numerical model was used to generate simulated observed LPs and nested sampling was used to infer spot locations using the analytical model.

For each trial, we reviewed the results in the form of corner plots and judged each trial as a pass, blend, or fail. A pass would be awarded if the inferred locations for both spots were within 1σ of the true values with 1 distinct peak; a blend would be awarded if the 1σ range for the two spots overlapped in either longitude or latitude; and a fail would awarded if either condition above was met. All cases passed the test for a grid with ΔL and ΔL from 0° to 60° with increments of 15°. For the case of an irregularly shaped spot, the analytical model would likely interpret this as either multiple overlapping individual inhomogeneities or one Gaussian spot with a reduced contrast (due to the portion of non-contrasted background being incorporated into the larger spot).

3.3. Photometric Validation

To truly replicate results in Karalidi et al. (2015), we needed to compare our photometric model to that shown in their work. To this end, we performed the same analysis as shown in Section 3.2 only with photometric data.

Again, we set the rotational speed to v sin i to 50 km s−1 and the viewing angle to nearly edge-on, i = 85°. The primary and secondary spots were at the same locations with identical spot radii and contrasts as outlined in Section 3.2. The numerical model was used to generate a photometric light curve as described in Section 2.3 while the analytical model was used to compute the likelihood within a nested sampling Bayesian inference. Both the analytical and numerical photometric light curves were smoothed using a Gaussian filter as described in Sections 2.2 and 2.3.

The results are very similar to those shown in Figure 5 with improved spot resolution versus that seen in Karalidi et al. (2015). Within the context of this sample case, the analytical code resolves both spots' locations for both spectroscopic and photometric data.

As in Section 3.2, we conducted an exhaustive exploration for varying longitude and latitude displacements for i = 30° and 60°. The spot size and contrast as well as the noise level (1%) were identical to those discussed in Section 3.2.2. Just as in that section, the photometric analytical model successfully retrieved every input scenario.

4. Application to Luhman 16B

With the analytical Doppler imaging method validated versus an numerical model, the next step is to perform an external validation with real-world observational data. Similar to Luger et al. (2021), we choose to validate our technique against Luhman 16B (WISE J104915.57-531906B).

4.1. Luhman 16AB System

Luhman 16AB are two of the solar system's nearest neighbors at 1.9955 ± 0.0004 pc, which also makes the two binary components the nearest BDs (Bedin et al. 2017). Luhman 16AB was originally discovered by Luhman (2013) using images from the Digitized Sky Survey, the Two Micron All-Sky Survey, and the Deep Near-Infrared Survey of the Southern Sky. Luhman 16AB is comprised of Luhman 16A, a ${34.2}_{-1.1}^{+1.3}\ {M}_{J}$ primary with an L7.5 spectral class, and Luhman 16B, a ${27.9}_{-1.1}^{+1.0}{M}_{J}$ secondary with a measured spectral class T0.5 (Burgasser et al. 2013; Kniazev et al. 2013; Garcia et al. 2017; Ammons & Garcia 2019).

Similar to many BDs at the L/T transition, Luhman 16B demonstrates high photometric and spectroscopic variability (Biller et al. 2013; Gillon et al. 2013; Burgasser et al. 2014; Buenzli et al. 2015a, 2015b). Although the primary, Luhman 16A, is also considered variable, it is to a smaller degree than the secondary, Luhman 16B (Biller et al. 2013; Burgasser et al. 2014; Crossfield et al. 2014; Buenzli et al. 2015a). Several studies over the last decade have investigated Luhman 16B's variability and periodicity. Gillon et al. (2013) deduced an 11% peak-to-peak near-infrared (0.9 μm) variability and 4.87 ± 0.01 hr period over a 12 night campaign. Burgasser et al. (2014) estimated a similar 13.5% near-infrared (1.25 μm) variability. Apai et al. (2021) used TESS to observe 100 rotations of the binary system, finding maximum peak-to-peak variability to be approximately 4% with longer period, 10% variability seen over 100 hr timelines. Apai et al. (2021) attributed the strongest periodicity signal (5.28 hr) along with a secondary periodicity (2.5 hr) to the combined effects of differential rotation and planetary waves in Luhman 16B's atmosphere. In line with the above results, Heinze et al. (2021) observed a 5 hr period and variability in the red continuum changing from 7% to 13% on consecutive nights, providing a relatively consistent measured period and high-variability.

4.2. Previous Global Maps of Luhman 16B

Global surface maps of Luhman 16B's surface inhomogeneities have been generated from both spectroscopic and photometric observations. Crossfield et al. (2014) inferred spot latitude, contrast, and size based on a 1-Spot model. They also created the first surface map of the BD by applying the maximum entropy method outlined in Vogt et al. (1987).

More recently, Luger et al. (2021) reprocessed the Crossfield et al. (2014) data and used the open-source Python code Starry (Luger et al. 2019) to generate a global surface map. The map created by Luger et al. (2021) is qualitatively similar to that of Crossfield et al. (2014) with both studies deducing an dark equatorial region and a bright polar spot.

Using photometric data from the HST, Karalidi et al. (2016) generated three and four spot surface maps for Luhman 16A and B. For Luhman 16B, a four spot model with temperature (contrast) as a free parameter was qualitatively different from the maps created by Crossfield et al. (2014) and Luger et al. (2021), but still contained near equatorial spots and brighter spots in the polar region (Karalidi et al. 2016). In the following subsections, we will use archival data from both of these studies to produce Luhman 16B surface maps with the analytical method described in this paper.

4.3. Preparing Spectroscopic Data

To infer a global map of Luhman 16B's surface, we follow a process similar to that used by Crossfield et al. (2014) with the same VLT/CRIRES data collected on 2013 May 5. Crossfield et al. (2014) initially obtained 56 spectra (2.288−2.345 μm) with 300 s exposure times and combined these into groups of four to create high S/N LPs via the LSD method. We obtained the Luhman 16B data via open online archival and used the reprocessed template and observed spectra from Luger et al. (2021) to similarly create LPs via the process outlined in Section 2.1. Referencing Section 2.1, the template spectrum created by Crossfield et al. (2014) and reprocessed by Luger et al. (2021) is used to create a line mask matrix (M). Equation (4) is then used to compute the LP for each epoch.

An important point to note is the retrieved latitude's sensitivity to factors such as LP element number (n) and the regularization parameter (Λ). We selected n = 93 LP elements to ensure a ΔRV ≈ 1−2 km s−1 while also sampling the RV range of ±3σ. Λ = 500 was qualitatively selected to attenuate noise in the retrieved LP while preserving the observed signal. Varying these two parameters can vary the retrieved latitude, but all resulting variations of the retrieved latitudes are consistent with those of Crossfield et al. (2014).

4.4. Spectroscopic Results

The resulting LPs can be seen in Figure 7 along with the line fit from our 1-Spot analytical model. Due to the retrieved spot's radius and relatively low contrast (compared to the validation examples), the affected region of the LP is broad and the subtle shifts in flux are difficult to detect by eye, necessitating a computational Bayesian method. Using nested sampling, we inferred spot longitude, latitude, radius, and contrast by comparing the LPs generated from observed spectra to our analytical model. As in Crossfield et al. (2014) and Luger et al. (2021), we set inclination to i = 70° and v sin i = 26.1 km s−1. We manually varied the linear limb darkening coefficient (epsilon) and the number of LP elements (n) to fit the observed data. We convolved the LP with a Gaussian profile to account for instrument and other sources of line broadening resulting in the modeled LPs in Figure 7.

Figure 7.

Figure 7. LPs for 14 epochs, representing one rotation of Luhman 16B. Dashed red lines represent observed spectra from Crossfield et al. (2014) while the solid blue lines denote the 1-Spot analytical model developed in this paper. Each epoch is accompanied by the appropriate time stamp in hours with the period being ≈5 hr. LPs are offset by 0.02 for comparison, but are normalized to unity within the code. The mid-phase LP deviations are difficult to detect visually due to the dominant spot's large radius and low contrast creating shallow flux variations, thereby, necessitating computational methods to retrieve. Note the analytical model has been convolved with a Gaussian instrument profile to account for the spectral resolution.

Standard image High-resolution image

For the 1-Spot model, dynamic nested sampling resulted in an equatorial dark spot (latitude = 29°) with radius = 37° and contrast = 0.10, which can be seen in Figure 8. These results can be compared to those of Crossfield et al. (2014) and Karalidi et al. (2016), and the photometric results in Table 1. Our results match well with the values given in Crossfield et al. (2014), which inferred a latitude ≤31°, radius = 33° ± 7°, and contrast = 0.12 ± 0.03. Longitude was not included in Crossfield et al. (2014) as it is dependent on the beginning epoch of observation.

Figure 8.

Figure 8. 1-Spot (spectroscopic) corner plot. The mean inferred values along with 1σ quantiles displayed on the top of each column. Contours denote the 0.5σ, 1σ, 1.5σ, and 2σ regions. Created with Dynesty (Speagle 2020).

Standard image High-resolution image

Table 1. Summary of the Model Comparison a

ModelModeLatitudeLongitudeRadiusContrast
Crossfield et al. (2014)Spectroscopic≤31°N/A33° ± 7°0.12 ± 0.03
1-SpotSpectroscopic $29^\circ \frac{+8.44^\circ }{-8.90^\circ }$ $-160.39^\circ \frac{+2.08^\circ }{-2.21^\circ }$ $36.79{^\circ }_{-6.73^\circ }^{+7.73^\circ }$ ${0.10}_{-0.03}^{+0.04}$
Karalidi et al. (2016)Photometric (J)−20° ± 12°113.47° ± 6.46°38.89° ± 0.8°−0.180 ± 0.013
 (fixed bright spots)45° ± 12°186.08° ± 3.65°24.45° ± 1.46°−0.180 ± 0.013
  72° ± 10°283° ± 13°35.54° ± 1.67°−0.180 ± 0.013
Karalidi et al. (2016)Photometric (H)−18° ± 8°107.27° ± 5.48°32.42° ± 0.8°−0.180 ± 0.013
 (fixed bright spots)51° ± 20°190.91° ± 4.83°25.02° ± 1.13°−0.180 ± 0.013
  74° ± 6°288.6° ± 9.4°36.30° ± 1.18°−0.180 ± 0.013
1-SpotPhotometric (J/H) $-23.53^\circ \frac{+15.72^\circ }{-9.89^\circ }$ $110.08{^\circ }_{-2.79^\circ }^{+3.30^\circ }$ $37.58{^\circ }_{-5.37^\circ }^{+4.86^\circ }$ ${0.25}_{-0.06}^{+0.07}$
2-SpotPhotometric (J/H) $-27.20{^\circ }_{-10.86^\circ }^{+18.22^\circ }$ $119.30{^\circ }_{-5.22^\circ }^{+5.92^\circ }$ $34.86{^\circ }_{-4.85^\circ }^{+3.59^\circ }$ ${0.28}_{-0.06}^{+0.05}$
   $47.19{^\circ }_{-32.52^\circ }^{+24.41^\circ }$ $\ast 231.68{^\circ }_{-20.44^\circ }^{+15.54^\circ }$ $28.20{^\circ }_{-7.44^\circ }^{+7.64^\circ }$ $-{0.14}_{-0.09}^{+0.06}$

a Note. ∗ Denotes where +360° has been added to negative longitudes to more easily compare to external results where the longitudinal phase goes from 0° to 360° vs. −180° to 180°

Download table as:  ASCIITypeset image

The 1-Spot model reproduces the MCMC results in Crossfield et al. (2014) and shows them to be a reasonable representation of Luhman 16B's dominant feature, the equatorial dark spot (Figure 9). For models with a higher number of spots, the spectroscopic analytical method did not converge on a likelihood, leading us to believe the 1-Spot model best describes the Crossfield et al. (2014) data.

Figure 9.

Figure 9. Luhman 16B LP deviations. (Left) Deviations from Crossfield et al. (2014) observed mean LP. (Middle) Modeled deviations for the inferred equatorial dark spot. (Right) Residuals computed by subtracting modeled LP deviations from observed LP deviations.

Standard image High-resolution image

4.5. Photometric Results

Luhman 16A and B photometric data were collected on 2013 November 8 by the HST/Wide-Field Camera with the G141 grism, which operates between 1.1 and 1.66 μm. The data was originally collected, processed, and reported by Buenzli et al. (2015b). Karalidi et al. (2016) used those observations, as well as those collected by Buenzli et al. (2015a) with the G102 grism (0.8−1.15 μm) on 2014 November 23, to construct three and four spot surface maps using the Aeolus code developed in Karalidi et al. (2015).

Here we will seek to infer surface spots from the Buenzli et al. (2015b) data using the analytical method described in Section 2.2. The J and H near-infrared bands from Buenzli et al. (2015b) result in very similar light curves, so for this paper these bands are averaged.

In Figure 10, the Buenzli et al. (2015b) J- and H-band data, along with our 1- and 2-Spot photometric models are plotted. The 2-Spot light-curve model more fully represents Luhman 16B's photometric behavior. Dynamic nested sampling results in a model with improved evidence (${\rm{\Delta }}\mathrm{ln}{ \mathcal Z }=13.8$) for the 2-Spot model over the 1-Spot model. This is strong evidence that the 2-Spot model better represents the data (Kass & Raftery 1995). A summary of the photometric results can be found in Table 1.

Figure 10.

Figure 10. Photometric light curves from Buenzli et al. (2015b) along with 1- and 2-Spot analytical models. The 1-Spot model includes an equatorial dark spot while the 2-Spot model includes both an equatorial dark spot as well as a midlatitude bright spot. Based on the results of Karalidi et al. (2016), the unmodeled trends in the light curve are potentially due to additional bright and dark spots combined with Luhman 16B's observed short-term timescale variability.

Standard image High-resolution image

Further increasing the number of spots in our model leads to lower comparative evidence (${\rm{\Delta }}\mathrm{ln}{ \mathcal Z }$), so we conclude the strongest evidence exists for the 2-Spot model. The equatorial dark spot and midlatitude bright spot inferred by the 2-Spot photometric model echo the results seen in the photometric results in Karalidi et al. (2016) as well as the Doppler imaging maps in Crossfield et al. (2014) and Luger et al. (2021).

A summary of both spectroscopic and photometric results, as well as comparisons to those of Crossfield et al. (2014) and Karalidi et al. (2016), can be seen in Table 1. Both the spectroscopic and photometric results reasonably match those obtained in Crossfield et al. (2014) and Karalidi et al. (2016), respectively. Even though the results from Karalidi et al. (2016) rely on fixed bright spots (negative contrast for this paper's nomenclature), our 1-Spot and 2-Spot analytical models still approximately infer the average spot location and size of the J and H bands in Karalidi et al. (2016). This matches our expectations as we averaged the J and H bands before performing Bayesian inference via nested sampling. Furthermore, our 2-Spot photometric model resulted in a Bayesian information criterion value of 41.4, matching the J-band 2-Spot model value achieved in Karalidi et al. (2016).

5. Discussion

5.1. Persistent Atmospheric Features

Multiple data sets at different wavelengths have explored the surface of Luhman 16B. From these studies, evidence has emerged that the BD possesses a persistent dark equatorial spot and potentially a bright mid-to-polar latitude spot.

Persistent atmospheric features for Luhman 16B can be seen in the results of this paper as well as those in Crossfield et al. (2014), Luger et al. (2021), and Karalidi et al. (2016). As discussed above, Crossfield et al. (2014) inferred a dark equatorial spot using MCMC. Using the maximum entropy method, Crossfield et al. (2014) also created a global map that included a dark equatorial spot as well as a bright polar feature. Luger et al. (2021) reprocessed the same spectroscopic data and inferred a similar map with those same features. Karalidi et al. (2016) identified the dark feature and labeled it the Possible Persistent Cloud Structure (PPCS-1), which may be similar to Jupiter's Great Red Spot. Not only did it appear in the HST J and H bands as well as a G102 grism light curve, it also appears in data presented in Gillon et al. (2013), which used the I + z TRAPPIST filter (Karalidi et al. 2016).

These features may be consistent with those in Apai et al. (2021), in which Luhman 16B's planetary waves, potentially in the form of zonal jets, were demonstrated over 100 hr of TESS observations. It is possible the persistent dark spot forms in such a zonal jet at Luhman 16B's equator while mid-to-polar spots form within higher latitude bands. Apai et al. (2021) also detected what they thought to be a long period signal from a notional polar vortex on Luhman 16A. The existence of such persistent features could help shed light on the overall atmospheric structure and dynamics of Luhman 16B and similar ultracool objects.

While there is evidence for persistent features, variability remains a recurring theme in observations of BDs. L/T transition variability is likely caused by multiple atmospheric layers and clouds with disparate opacities. The variability of Luhman 16B (Gillon et al. 2013; Crossfield et al. 2014; Buenzli et al. 2015a, 2015b) may also arise from the different observational techniques and wavelengths considered. Not only did these studies observe Luhman 16B at different times, but they also probed different pressure levels and chemical species. The rich inhomogeneous data sets can complicate the picture when attempting to quantify the true variability of the object. However, they also offer an opportunity to understand the BD's three-dimensional atmospheric dynamics.

5.2. Doppler Imaging as Model Validation

Over the past decade, there has been extensive work in modeling atmospheric dynamics in lightly irradiated BD atmospheres through general circulation models (GCMs). Several drivers for robust circulation in the stratified overlying atmosphere of these BDs have been explored in the literature, including convection driving the stratified upper layer and cloud-radiative feedback. Showman et al. (2020) and Zhang (2020) both include comprehensive reviews of these efforts.

Convection from the interior perturbing the overlying stratified atmospheric layer has been a mechanism identified as a potential primary driver of BD circulation (Freytag et al. 2010; Showman & Kaspi 2013). The strength of these perturbations, along with the object's radiative timescale and drag at the radiative-convective boundary, appear to control the latitudinal extent and strength of jet stream formation (Zhang 2014; Showman et al. 2019; Tan 2022).

Cloud-radiative feedback has also been demonstrated as a plausible mechanism for driving atmospheric circulation. One-dimensional feedback between clouds, radiation, and atmospheric mixing is explored in Tan (2019). More recently, modelers have focused on three-dimensional feedback models in which cloud-radiative dynamics induce active atmospheric circulation thereby perpetuating cloud patchiness (Tan 2021a, 2021b; Lefèvre et al. 2022).

Placing Luhman 16B's large equatorial dark spot and midlatitude/polar bright spots within the context of a model such as that discussed in Zhang (2014) offers a means to connect Doppler imaging to atmospheric dynamics. In that work, the atmospheric structure is essentially a product of forcing from the convective interior and the radiative timescale of the overlying stratified atmosphere. Strong convective forcing and weak radiation lead to large equatorial jets with midlatitude to polar vortices and turbulence. It is possible that the equatorial dark spots are persistent clouds or storm systems within a jet stream band. The bright spots may be regions of lessened opacity and cloud thickness created by turbulent features in which we can observe deeper layers of the BD's atmosphere.

Combining BD atmospheric circulation models with Doppler imaging techniques could potentially lead to more accurate surface maps as well as validation for the development and selection of GCMs. Follow-on observations of Luhman 16B and other ultracool targets may result in more accurate atmospheric models.

5.3. Doppler Imaging and Chemical Disequilibrium

Doppler imaging focused on specific chemical species' spectral lines could help to further illuminate the mechanisms at play in BD and gas giant exoplanet atmospheres. Tracing specific species such as CO, CH4, NH3, and molecular nitrogen (N2), provides insight into atmospheric circulation. At temperatures above the L/T transition, carbon's chemical equilibrium state transitions from favoring CO (above 1470 K) to CH4 (below 1470 K) (Fegley & Lodders 1996). A similar chemical equilibrium transition occurs for nitrogen. Below 630 K, NH3 should be the dominant species while N2 is favored by chemical equilibrium above this value (Fegley & Lodders 1996).

Chemical disequilibrium in ultracool atmospheres provides evidence for robust circulation. The expectation for BD spectra at temperatures at or below the L/T transition is to favor higher relative abundances of CH4 and NH3 and lower abundances for CO and N2. However, numerous BD atmospheric observations dating back to Gliese 229B have detected above equilibrium values of CO (Noll et al. 1997; Oppenheimer et al. 1998) and a dearth of NH3 (Saumon et al. 2000, 2006).

This disequilibrium was expected by atmospheric modelers based on observations of Jupiter and Saturn (Fegley & Lodders 1996). Chemical disequilibrium might indicate atmospheric circulation with convection or advection transporting CO into the higher stratosphere while sequestering NH3 within the BD and gas giant interiors (Noll et al. 1997; Griffith & Yelle 1999; Saumon et al. 2000; Leggett et al. 2007).

Applying Doppler imaging to CO, CH4, and NH3 spectral lines could potentially create BD and gas giant exoplanet surface maps detailing regions of increased convective circulation as well as other types of vertical transport observed in solar system stratospheres. Connecting these maps to dynamical models would allow for a greater understanding of ultracool objects' atmospheric structure and dynamics.

6. Summary

We have developed an analytical framework for inferring surface inhomogeneities on ultracool objects. The analytical framework is internally validated versus a numerical model for different inclinations and varying latitudinal and longitudinal displacements in both spectroscopic and photometric modes. This model which employs nested sampling demonstrates a computationally inexpensive and easy-to-implement approach to inferring surface inhomogeneities and is capable of retrieving spot location, size, and contrast with relatively high resolution. We then applied the analytical technique to archival spectroscopic and photometric data and retrieved similar results to those found in the literature. Moving forward, this technique can be useful for Doppler imaging the surfaces of increasingly faint ultracool objects and those with spectra made available by upcoming ELTs. Additionally, this technique could be paired with GCMs and theoretical models to discover and validate the complex dynamics driving atmospheric circulation in ultracool objects.

The authors would like to thank the United States Air Force Academy, Department of Physics, for sponsoring the graduate work of the first author. Special thanks to the anonymous reviewer for providing insightful and constructive feedback. We would also like to thank Dr. Rodrigo Luger of the Flatiron Institute Center for Computational Astrophysics for advice on modeling stellar surfaces. Additionally, we would like to thank Dr. Theodora Karalidi for sharing the Buenzli et al. (2015b) photometric data used in this article. Finally, we would like to thank the Group for Studies of Exoplanets (GFORSE) at The Ohio State University, Department of Astronomy, for continuous feedback throughout the development of this research.

The views expressed in this article are those of the author and do not necessarily reflect the official policy or position of the Air Force, the Department of Defense, or the U.S. Government.

Software: Dynesty (Speagle 2020), SciPy (Virtanen et al. 2020).

Please wait… references are loading.
10.3847/1538-4357/ac75b9